function lambda_difference_pt = underinvestment(oParams,bParams,setts,f,ptsim,ctsim,lambdasim)
% Underinvestment
% This function was created to calculate the level of underinvestment
% relative to first-best

% Variables:
% bParams,oParams, setts	see readme
% f                 rate of creative destruction
% ptsim      		simulated values for the number of product lines p
% ctsim      		simulated values for the coupon p
% lambdasim      	simulated values for the investment intensity lambda

% Set parameters
	pbar = oParams.pbar;
    pi = oParams.pi;
	xi = bParams.xi;
    eta = bParams.eta;
    
% Calcualte UpMat
    upMat = updatematrix(bParams,oParams,setts);

% Calculate first best equity value
    c_vec = unique(ctsim);
    n_c = length(c_vec);
    lambda_first_best = zeros(n_c,pbar+1);
    for i = 1:n_c
        setts2 = setts;
        setts2.C = -((pi-xi)/(1-pi)) * c_vec(i);
        if (pi-xi)* c_vec(i) > (1-pi) * eta
            disp('There is no default at zero and therefore a problem with the code');
        end
        [~, ~, lambda_first_best(i,:), ~, ~, ~, ~] = equity_value_coupon(bParams,oParams,setts2,upMat,f);
    end
    
% Determine difference in innovation intensities
    lambda_difference_sim = ptsim;
    for p = 0:pbar
        for i = 1:n_c
            lambda_difference_sim(ptsim==p & ctsim==c_vec(i)) = (lambda_first_best(i,p+1) - lambdasim((ptsim==p & ctsim==c_vec(i))))./lambdasim(ptsim==p & ctsim==c_vec(i));
        end
    end
    
% Average over c (conditional on p)
    lambda_difference_pt = zeros(1,pbar);
    for p = 1:pbar
        lambda_difference_pt(p) = mean(lambda_difference_sim(ptsim==p));
    end
end